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ABSTRACT 


This  paper  shows  that  partial  differential  equations  may 
be  a  possible  area  of  application  of  mathematical  programming 
The  solution  of  Lapiace's  equation  with  Neumann's  condition  is 
shown  to  be  a  minimum  cost  network  flow  problem  with  cost  pr 
portional  to  the  arc  flow.  An  algorithm  of  solving  minimum 
quadratic  cost  network  flow  is  given. 


I.  Introduction: 

The  purpose  of  this  paper  is  to  present  a  new  numerical 
method  for  solving  Laiiace's  equation  with.  Neumann's  condition 
prescribed  cn  the  boundary.  This  method  uses  the  idea  in  the 
theory  of  network  t  ews  whicn  is  a  branch  cf  mathematical  program¬ 
ming  usual!/  considered  to  1;  unr  elated  to  the  subject  of  partial 
differential  equations-  For  6 ' mplic ity  of  exposition  we  shah,  dis¬ 
cuss  the  case  in  two  dimensions  The  generalization  to  three  dimen¬ 
sions  is  obvious 

Consider  a  simply  •  c  nnccted  dorrciin  G  in  the  plane  whose 
boundary  T  is  a  simple  closed  curve  We  are  interested  in  find 
ing  the  solution  of 
(1)  V  0  in  G 


0d>  3d) 

with  tt-I-  as  prescribed  or  1  ,  where  .  , 

on  r  on  is  the  normal  deriva¬ 

tive  of  <f>  consider  =d  to  be  positive  if  <f>  decreases  its  value 
across  the  boundary  from  the  outside  regicr  to  C-  In  order  to 

have  a  unique  so.uticn,  r  is  necessary  and  sufficient  to  prescribe 
9<*> 


on  T  such  that  the  ,'.:.e  Integra: 


(2) 


This  prob.em  arises  :  -  om  many  problems  of  mathematical 
physics;  we  shau  only  give  the  physical  interpretation  from  the  point 
of  view  of  fluid  mechanics  We  consider  this  problem  as  that  of  in¬ 
compressible  ir  rotations1  flow  and  thus  is  the  velocity  of  the 

9$  _ 

flow  where  -  on  i  l -s  ini  lux  or  outflux  cf  the  liquid  across  the 


boundary  where  (i)  .s  the  equation  of  continuity  Since  there  is  no 


source  or  sirk  inside  C  ,  \Z)  expresses  the  condition  that  the 
total  inflow  must  equal  the  total  outflow.  We  shall  call  the  parts  of 

o  u 

F  for  which  are  positive  ime  sources,  and  the  part  of  T 

for  which  -jj—  are  negative  line  sinks. 

The  usual  numerical  method  is  to  replace  differential  equation 
(1)  by  difference  equation  with  uniform  or  non-uniform  grids,  and 
then  calculate  the  value  o'.  the  function  at  the  discrete  points  where 
the  spacing  between  grid  uies  is  determined  by  the  degree  of  accu¬ 
racy  de sir  ed. 

For  a  non-uniform  rectangular  grid  (see  Fig.l)  ,  the  differ¬ 
ence  equat’cn  icz  the  Laplace's  equation  at  point  P  becomes 


hE^E  '  K\V; 


9P  "  *P 

V  ‘  £,  vV'  N  N  S' 


*P  "  *s 


S'  N  S' 


F  ie  • 


l 


4>e  ■  4>p 

Denoting  ^ ^ — j-  ,  etc  ,  we  can  rewrite  (3)  as 

EE  W 


(4) 


fEP  '  fPW  +  fNP  "  fPS  "  °  ‘ 


If  we  consider  the  equation  at  the  point  S  ,  using  this  notation,  we 
also  have 


,c,  ,  't’p  "  ‘•’s 

‘  ps  ’■  hs(hs  +  hT)  - 

Since  h^  h^,  in  general,  the  symbol  fpg  etc.  ,  seems  to  be 
not  uniquely  defined.  We  shall  first  assume  that  all  h's  are  equal 
and  then  discuss  the  non-uniform  grid  case  later.  For  the  case,  all 
h's  are  equal,  and  the  f.  are  just  the  differences  between  func¬ 
tions  at  points  divided  by  2h^  and  are  uniquely  defined.  Equation 
(4)  suggests  a  strong  analogy  between  the  difference  equation  and  the 
physical  model  of  water  flowing  in  pipes. 

We  shall  now  consider  a  physical  model  for  the  grid  system  to 
be  a  network.  The  boundary  value  problem  of  a  network  has  already 
been  discussed  in  Duffin  [  8]  and  Birckhoff  and  Diaz  [  2]  .  As  usually 
done  in  the  theory  of  network  flows  (see  the  most  complete  descrip¬ 
tion  by  Ford  and  Fulkerson  [  9]  )  we  consider  a  network  which  con¬ 
sists  of  nodes  N.  and  arcs  B..  connecting  N.  and  N.  . 

i  ij  i  J 

The  value  of  the  flow  from  N.  to  N.  along  the  arc  B..  is  de- 

i  J  ij 

noted  by  f. .  and  we  have  f. .  s  -f .. 

ij  V)  P 

The  conservation  of  liquid  at  a  node  N.  is  expressed  by  the 
equation 

(6)  V  f . .  -  Vf.  .  =  0  for  all  i  £  s,t, 

L  u  L  kl 

j  k 
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where  f  >0  and  f.  .  >  0 
’J  —  ki  — 

If  a  node  N  is  a  source  node  of  strength  v  ,  we  have 
s 

<7>  I£SJ  =  * ; 

j 

and  for  a  sink  node  N  of  strength  v  ,  we  have 


1 


We  have  now  replaced  a  continuous  domain  by  a  discrete  net¬ 
work  with  line  sources  on  F  replaced  by  source  nodes,  and  line 
sinks  by  s,nk  nodes.  The  strength  of  a  source  representing  a  short 
line  source  y  is  determined  by 

(9)  V  =  j  |4dS 

JY 

and  similarly  for  a  short  line  sink.  In  the  following  discussion,  we 
shall  assume  that  there  is  only  one  source  node  and  one  sink  node. 

The  reader  familiar  with  network  flow  problems  or  problems 
of  operations  research  will  suspect  that  this  problem  becomes  a 
minimum  cost  flow  problem  [  9]  or  that  of  Hitchcock  transportation 
problem  [  1Z ]  with  the  sources  as  supplies  and  sinks  as  demands. 

But  this  is  not  a  standard  minimum  cost  flow  problem,  as  shown 
later-  For  one  thing,  the  network  with  infinity  branch  capacities 
can  have  many  minimum  cost  flow  patterns  in  the  network  if  the  cost 
is  defined  to  be  proportional  to  distance,  where  the  solution  of 
Laplace's  equation  is  known  to  be  unique. 
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II  Network  Formulation: 


To  solve  this  question,  we  now  turn  to  Dirchlet1  s  principle 
(see  the  book  by  Courant  [5])  .  The  Dirchlet  integral  of  a  function 
g  defined  on  C  is  given  by 


,  2  2 

'gx  +  V  dX  dy 


where  g  denotes  the  partial  derivative  of  g  with  respect  to 
x  ,  similarly  fcr  g^ 

Now  we  aucte  from  Courant  [5]  . 

Dirchle^s  principle:  Given  a  domain  G  whose  boundary  T  con¬ 
sists  of  Jordan  curves.  Let  g  be  a  function  continuous  in  G  +  T  , 
piecewise  smooth  in  G  ,  and  with  finite  Dirchlet  integral  D[g] 
Consider  the  crass  cf  a*;  functions  <j>  continuous  in  G  +  T  , 
piecewise  smooth  in  G  ,  and  having  the  same  boundary  values  as 
g  .  Then  the  problem  cf  finding  a  function  4>  for  whic  h  D[4>] 
attains  a  minimum  d  ha  s  a  unique  solution  4>  =  u  This  func 
tion  u  is  the  solution  of  the  boundary  value  problem  fcr  V^u  -  C 
with  the  prescribed  boundary  vaiue  g  on  T 

Dir  chief's  principle  is  certainly  correct  if  in  the  boundary 


value  T.'rob’wm,  the  normal  derivative  of  g  is  prescribed  instead 
of  e  The  finite  ana'ogy  of  (10)  for  a  discrete  network  is  then 


(11) 


D[4>] 


N.  P. 


(The  classical  approach  is  to  differentiate  (11'  with  respect  to 


and  solve  the  simultaneous  linear  equations  thus  obtained  ) 


5 


For  the  case  all  h  -  1 


we  have 


(IZ)  DM  =  Dt'ij]  =  j  £  ffj 

i.  j 

So  the  network  problem  becomes  that  of  minimizing  (1Z)  sub¬ 
jected  to  the  condition  (6)  (7)  (8)  .  This  is  a  mathematical  program¬ 
ming  problem  with  linear  constraints  and  non-linear  objective 
function  (1Z)  .  With  the  solution  space  being  convex  and  the  objec¬ 
tive  function  quadratic  (see  for  example  [lb]),  relative  minimum  in 
the  solution  spare  implies  absolute  minimum.  Furthermore,  the 
solution  is  unique. 

Instead  of  solving  the  problem  as  a  quadratic  programming 
one,  we  shall  develop  a  network  flow  method.  For  convenience,  we 
assume  that  a  scale  system  is  chosen  such  that  the  strength  of  Ng 
and  of  are  very  large  integers  and  "one"  is  the  smallest  posi¬ 

tive  constant 

Let  us  define  the  cost  c^  of  shipping  one  unit  of  flow  from 

N.  to  N  as 
i  J 


(13) 


c . . 

ij 


r  f.  .  if  f.  .  >  0 

U  - 

-f..  if  f.  >  0  . 

Jl  Jl 


This  problem  is  now  a  minimum  cost  flow  problem  with  the  cost  of 
shipping  along  an  arc  depending  on  the  direction  and  the  amount  of 
flow  in  that  arc 

We  are  interested  in  the  flow  pattern  which  satisfies  the  sup¬ 
ply  and  demand  conditions  as  defined  by  (6)  (7)  and  (8)  with  the  total 
cost 
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(14) 


D 


f.. 


'J 


a  minimum.  Let  us  call  this  minimum  cost  flow  pattern  the 
optimum  flow  pattern. 

Clearly,  in  this  optimum  flow  pattern,  if  a  set  of  arcs 

B.,,B,-, . B  form  a  cycle,  then  we  cannot  have 

ll  12  mi  1 

f.,,f,,  all  positive.  Otherwise  we  can  subtract  from 

ll  12  mi 

each  arc  the  amount  <5  =  min  (f. ,,  f f  .)  >  0  ,  with  (6) 

i  l  1 2  mi 

(7)  and  (8)  still  satisfied,  but  with  the  value  of  (J4)  reduced.  This 

is  no  surprise  as  we  know  from  the  continuous  case,  the  solution 

of  Laplace's  equation  represents  irrotational  flow.  This  will  be 

called  the  condition  of  no  circulation. 

Consider  now  the  case  for  which  there  is  a  cycle  formed  by 

arcs  in  which  we  pick  a  node  N.  and  a  node  N.  .  Let  us  con- 

i  J 

sider  the  flows  f . f f  in  the  first  path  and  the  flows 

il  12  mj  r 

f.  ii  ifii  f  , .  in  the  second  path.  Note  these  arc  flows 

1 1 '  1' ,  2 '  m'j  r 

may  be  positive  or  negative.  Then  we  must  have 
<i5>  = 


where  the  left  side  represents  the  marginal  cost  of  shipping  one 
unit  of  flow  along  the  first  path,  and  the  right-hand  side  represents 
the  marginal  cost  of  shipping  along  the  second  path.  If  (15)  is  not 
true,  we  can  ship  one  unit  along  the  cheaper  path  and  go  back  on  the 
more  expensive  path.  This  will  make  (6)  (7)  and  (8)  all  satisfied 
but  make  (14)  reduced. 
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We  can  now  state  the  algorithm  of  constructing  the  optimum 

flow  pattern  This  algorithm  is  similar  to  the  algorithm  used  in 

[  3]  ,  [9]  and  [  1 5 ]  f  or  constructing  minimum  cost  flow  where  the 

cost  is  the  same  as  the  length  of  the  arc  or  as  its  negative 

Let  us  define  the  cost  of  a  path  from  N  to  N  as  the  sum 

s  t 

of  c o 8 1 s  of  arcs  travelling  along  that  direction. 
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Ill 


The  Method  and  its  Proof: 


This  algorithm  can  be  described  as  follows.  Assume  that  we 

have  sent  i  units  of  flow  from  N  to  N  (i  *  0  1,  .  ,  ,  v!  . 

s  t 

Step  I:  Find  the  cheapest  path  from  Ng  to  N  with  the 

costs  cf  arcs  defined  by  (13',  Let  the  path  cost  be  C.  (i.  e  .  the 

i  +  i 

t  H 

cost  of  the  (:+!’  ‘  path). 

Step  2:  If  C.  •-  C  ,  ^  0  ,  go  to  Step  3.  If  C.  -  C.  ,  =  0 

or  C.  <  C.  ,  ,  send  one  unit  of  flow  from  N  to  N  a.onp 

1  l+l  s  t  B 

the  cheapest  path.  Redefine  the  costs  of  arcs  on  the  path.  Go  back 
to  Step  1. 

Step  3:  Let  f.^  be  the  arc  flow  from  ISL  to  ,  then 

4>.  -  <t>,  =  . 

J  k  jk  i 

There  exist  many  papers  about  the  algorithm  of  finding  the 
cheapest  path  (or  the  shortest  path).  We  shall  only  mention  the 
first  paper  on  this  subject  by  Ford  [9]  ,  the  later  modification  by 
Dantzig  [7]  ,  the  matrix  formulation  by  Berge  [l]  ,  and  the  case  of 
many  sources  and  many  sinks  to  the  technique  by  Gomory  and  Hu 

[Hi 


We  shall  sav  N  covers  N.  if  there  is  a  positive  flow  f.  . 

1  J  ij 

from  N.  to  N.  .  Since  there  is  no  circulation  in  the  flow  pat- 
i 

tern,  we  can  consider  the  set  of  all  nodes  to  be  partially  ordered. 
For  the  network  with  its  nodes  partially  ordered  by  its  flow  pattern 
we  can  also  find  the  most  expensive  path  (without  cycles)  from  Ng 
to  N  (see  to*  example  [9l).  For  the  optimum  flow  pattern,  the 
costs  of  the  most  expensive  path  and  the  cheapest  path  should  be 
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equal 

First  we  prove  that  this  algorithm  does  give  a  flow  pattern  in 
which  costs  of  ail  flow  paths  are  equal.  Assume  that  there  are  two 
paths  and  with  the  cost  C  (P^)  <  C  (P?)  .  In  sending 

the  last  unit  flow  along  P^  .  we  are  not  using  the  cheapest  path 
contradictory  to  the  assumption  that  we  always  use  the  cheapest  path. 

Second,  we  prove  that  when  costs  of  all  paths  are  equal,  the 
optimum  solution  is  obtained.  Assume  that  there  is  an  optimum 
solution  with  the  flow  pattern  different  from  that  obtained  by  the  algo¬ 
rithm.  Assume  the  path  cost  of  the  pattern  obtained  by  the  algorithm 
is  C(h)  and  that  of  the  optimum  solution  is  C(C)  . 

Since  it  is  necessary  for  an  optimum  solution  to  have  all  paths 
the  same  cost,  if  the  optimum  solution  is  not  obtained  by  this 
algorithm,  it  must  have  all  paths  the  same  cost  with  C  (0)  <  C  (h)  . 
Consider  in  the  stage  of  the  algorithm  for  which  we  have  shipped 
v  -  1  units  of  flows  Then  we  could  send  the  last  unit  along  a  path 
with  cost  C  (0}  and  reduce  the  total  cost  Since  this,  by  assump¬ 

tion,  is  not  possible,  this  means  the  Zfj^  in  the  path  is  greater 
than  Z  f ?.  in  that  Dath  and  also  for  any  part  of  the  path.  This 

hi  G 

means  f.  .  >  f  If  this  is  true  for  all  paths,  since  optimum  so- 

U  U 

lution  gives  v  units,  then  by  (7)  this  means  we  have  shipped  more 
than  v  units  of  flow,  contradiction. 

In  order  to  speed  up  the  convergence  of  this  algorithm,  lots  of 
rules  can  be  given  for  sending  flows  from  the  source  to  the  sink 
when  there  are  many  arcs  with  f . .  -  0  in  the  arcs.  We  shall  not 
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discuss  such  rules  here. 


The  case  of  many  sources  of  strength  v.  and  many  sinks  oi 
strength  can  be  taken  care  of  as  usually  done  in  network  flows 

Also,  an  upper  bound  on  the  number  of  steps  to  solve  this  problem 
can  be  given  ( see  f 15 ] ) . 
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IV.  Examples: 


f 


11,  15 


Consider  Fig.  2  where  it  is  assumed  that 
=  8  and  ^  =  0  everywhere  else,  and 
According  to  the  algorithm,  we  first  send 


=  8 


h  =  1 


(As 

C3,  7 


as  c 


B 


3.7 


c 


7.11 


Step  1  The  cheapest  path  from  N^  to  costs  zero. 

f3. 7  "  f7,  ii  =  fn,  15  ~  0  imtially- ) 

Step  2:  Send  one  unit  of  flow  along  the  path.  Now 

^  c7.11  =  cll,  15  =  1  and  C 7 ,  3  =  C 11 ,  7  =  C 1 5 , 11  =  -1  • 

Step  h  The  cheapest  path  from  to  costs  now  2 

3,  7  =  1  ’  c6.  7  =  °  ’  C6.  10  =  0  ’  c10, 11  =  °  '  CU.  15  =  1  • 

Step  2.  Cq  =  0  <  Cj  =  2  .  Send  one  unit  of  flow  along  the  arc 


*  B7,  6  ’  B6,  10 

o 

CQ 

and  BU,  15 

Now 

=  1  ,  cn  ]5  =  2  , 

c7,  6  “  1 

'  c6, 10  =  1  • 

C10,  11 

By  this  algorithm,  we  will  send  two  more  units  along  the  path 


B3,  7  *  B7,  11  ’  Bll,  15  Wlth  C3,  7  ~  4  ’  C7,  11  _  3  *  cll,  15  “  4 


C7,  6  ~  1  1  c6,  10  "  1  ’  C10,  11  “  1  Wlth  C4  =  8 


Then  =  8  =  £  0  .  We  would  go  to  Step  3  and  obtain 


the  differences  between  functions. 

For  example,  4>-,  -  =  2  •  l3  •  3  •  (  ^ )  =  12 

The  {..  are  shown  in  Fig.  3  . 
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Lt*  as  ooriSidei  another  txampie.  A  special  example  such 
that  tw  very  Simple  consideration,  no  iteration  or  solution  of 
simultaneous  equations  is  recessary. 


Consider  Fig.  4  where  we  assume  that  all  normal  derivatives 

9d)  9(f> 

are  zero  except  that  we  have  -  14  at  Nj  and  =  -14  at 

N  ,  Due  tc  the  symmetry  of  this  configuration,  we  can  only  con¬ 
sider  the  pcints  N  N )  and 

Assume  the  amount  of  flow  from  N  ,  to  Nz  is  x 

4  6 

Then  from  the  :  onside:  a'ior.  of  flow  at  N,  we  have  i.,  -  f ,  . ,  -  x. 

b  46  64' 


Also,  f^-  and  f-^.  are  ecua ;  by  symmetry  Furthermore, 

f ,  r-  +  fc,  f.  /  +  i,  a  ■  otherwise  a  reduction  in  cost  can  be  a- 
45  54  46  64 

chieved.  Therefore  f,,  _  ic  A  ,  -x  Clearly  f , .  -  f . ,  +  f  ,  _  =  2x  , 

45  54*  7  14  46  45 

and  from  symmetry  and  the  equal  path  cost  condition 

hi  =  lii  B  4(<24  ♦  '45*  2  ~Z  •  ™s  will  make  =f23  +  f24  = 


3x  ,  7x 

T  +  Zx  "  T 


With  4j-T  =  14  =  I  +  f  2.  =  7 x  We  let  x  =  2 


and  get  the  differences  between  all  p»oints  as  shown  in  Fig.  5. 


.4 


V.  Discussion 


A  few  comments  are  in  order. 

1  This  paper  treats  the  Laplace's  equation  from  the  point  of 
view  of  mathematic  a.l  programming.  The  potentials  at  nodes  are 
dual  variab.t.s  tc  the  ws  in  arcs,  and  th^y  play  the  same  r  cl  e  as 
prices  in  the  •ar.gr3k’fc  o*  linear  programming 

2.  For  Poisson  s  equation  -  v(x,  y<  ,  the  same  inter¬ 

pretation  of  irrctationa!  incompressible  flow  with  sources  and  sinks 
on  the  boundary  and  inside  the  region  G  can  be  given.  The 
Dirchlet  integral  X  should  be  replaced  by 


-f 


2v  .‘x,  y.  4>)  dx  dy 


And  the  ccnditior.  2 


)  on  the  boundary  should  be  replaced  by 


d.s  + 


v ; x, y  dx  dy 


C 


Although  i  variational  formulation  of  a  general  elliptic 

partial  differential  equation  of  second  order  is  possib'e.  whether  the 

method  presented  here  is  applicable  is  unknown 

1  For  Dir  h.et  s  condition  prescribed  cn  the  boundary,  we 

have  just  the  same  number  of  unknowns  as  the  nimbei  of  equations 

The  same  is  'r^e  lor  Neumann's  condition  it  we  iix  4>  at  one  point. 

If  we  introduce  f.  as  variables,  then  we  double  the  number  of 

U 

unknowns  The  solution  will  not  be  unique  unless  we  use  the 
Dirchlet' s  Integra,  as  the  objective  function 

4  The  netwoiK  used  in  the  theory  of  network  flows  is  quite 
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arbitrary  Here  the  network  is  quite  regular  and  there  is  no 
capacity  constraints  on  the  arcs.  A  better  numerical  method  can  be 
developed  specialty  for  this  kind  of  network  Also  the  unimodular 
property  ,  see  for  example  [4]  and  [14]  ;  of  the  incidence  matrix 
describing  the  r.etwcik  gives  certain  insight  to  the  network  flow 
problem  This  may  plav  an  important  rote  in  further  development 
a‘ong  the  une  of  rht  numerical  methods. 

5  When  the  grids  are  not  uniform,  we  have 


This  makes  the  flow  in  an  arc  not  unique  But  this  can  be  taken  care 


ol  by  introducing  a  mu  tipiier  on  the  arc  from  Np  to  Ng  . 
This  muJLtip.irr  depends  only  on  the  geometrical  configuration  of 
the  grid.  The  paper  t  y  Jeweii  [14]  on  networks  with  gains  can  be 


applied  Fight  now,  few  papers  have  been  written  on  networks  with 


gains . 


-  IB 
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